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■ We analyzed the stochastic behavior of systems controlled by autocatalytic reaction A + X ^ 

' X + X, X B provided that the distribution of reacting particles in the system volume is 

I uniform, i.e. the point model of reaction kinetics introduced in arXiv:cond-mat/0404402 can be 

^SJ . applied. Assuming the number of substrate particles A to be kept constant by a suitable reservoir, 

we derived the forward Kolmogorov equation for the probability of finding n — 0,1, ■ ■ ■ autocatalytic 
particles X in the system at a given time moment. We have shown that the stochastic model results 
in an equation for the expectation value of autocatalytic particles X which differs strongly from the 
kinetic rate equation. It has been found that not only the law of the mass action is violated but 
' also the bifurcation point is disappeared in the well-known diagram of X particle- vs. A particle- 

, concentration. Therefore, speculations about the role of autocatalytic reactions in processes of the 

"natural selection" can be hardly supported. 
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Introduction 



^ ■ As in the previous paper we define the system S as an aggregation of particles capable for autocatalytic reactions. 

[ Symbols X and A are used for notations of autocatalytic and substrate particles, respectively. We denote by B the 

' , ' particles of end-product which do not take part in the reactions. We assume that the system is "open" for particles 

^ \ A, i.e. the number of particles A is kept constant by a suitable reservoir. However, the system is strictly closed for 

Q . the autocatalytic particles X. 

^ ' In this paper we investigate systems which are governed by reactions 

A + X^X + X, X + X^A + X, X^B, (1) 

' where kA, kx and kj are the rate constants. These reactions became interesting since convincing speculations were 
. published about the "natural selection" based on their properties. 

The organization of the paper as follows. After brief discussion of the kinetic rate equation, in Section I we derive 
and formally solve the forward Kolmogorov equation for the probability p„(i) of finding n = 0, 1, . . . X particles in 
• a system of volume V at time moment t > 0. Defining the conditions of the stationarity we determine the stationary 
I probability limt^oo Pn{t) = Wn, and analyze its properties. In Section II we present a modified stochastic model of 
, the autocatalytic reactions ^ which is capable of reproducing the solution of the rate equation, however, brings 
' about such a large fluctuation in the stationary number of X particles, that the mean value loses practically all of its 
^ , information content. In Section III we define the lifetime of a system controlled by reaction and calculate exactly 
^ • the extinction probability, as well as the mean value of the system lifetime. Finally, in Section IV we summarize the 
main conclusions. 



I. STOCHASTIC MODEL 

In order to make comparison, first the well-known results of the kinetic rate equation are briefly revisited, and then 
the stochastic model of the reactions ^ will be thoroughly analyzed. 

A. Rate equation 

Let m{t) be the number of X particles in the volume V of the system S at the time moment t > 0, and denote by 
(t) = m{t)/V the number-density of the X particles and by ca that of A particles which is kept strongly constant 



c 



'Electronic address |lpal@rmki . ktki . hu| 



2 



by suitable reservoir. According to the kinetic law of mass action we can write 

dcit) 



dt 

where 



= kACA c{t) - kx c'{t) - kd c{t) = kx c(t)[cR - c(t)], (2) 



kACA - kd 

CR = — (3) 



is the critical parameter of the reaction. Taking into account the initial condition c(0) = cx we obtain immediately 
the solution of ^ in the form: 

,,^ CxCRe^Y>{cRkx t} / n j U\ -f O IA\ 

c{t) = -f 7 — 7— -r TT, if CR ^ 0, and c{t) = — — -, if cr = 0. (4) 

CR + cx(exp{cRkx t} ~ 1) 1 + cxkxt 

Clearly we have two stationary solutions, namely 

ciJ' = CR, if CR > and c^^J = 0, if cr < 0. (5) 

It is easy to show that the solution of c^l'' is stable, if ca > ka/kA — Cbf, i.e. if cr > 0, while cff' is stable, if 
ca < kd/kA = Cbf, i.e. if cr < 0. Let Sc a small disturbance in c. It follows from Eq. ||2Jl that 

dSc 

~ kx [cr - 2c(t)] Sc, 

and so, we see immediately, if c{t) c^J' = cr, then Sc « exp{~kxCRt}, i.e. c^P is stable, when cr > 0. Similarly, if 

(2) (2) 

c{t) clf — 0, then 5c « exp{kxCRt}, i.e. c^j is stable, when cr < 0. 

The stationary density of X particles versus ca can be seen in FIG. 1. It is clear that ca = ctf is a "bifurcation 

point", since if ca > Cbf then there are two stationary solutions but among them only one, namely the solution c^^^ 
is stable. By decreasing the density ca adiabatically below the critical value Cbf, the autocatalytic particles X are 
dying out completely, and it is impossible to start again the process by increasing the density ca above the critical 

Cbf. 



Cat 




Cbf 

FIG. 1: Dependence of the stationary number-density of X particles on the number-density c^. The thick lines refer to stable stationary 
values. 

In the next subsection we will analyze the stochastic model of reversible reactions JQ). It will be shown that in the 
stochastic model the long time behavior of the number of autocatalytic particles is completely different from that we 
obtained by using the rate equation 

It is to mention that interesting and seemingly convincing speculations were published about the possibility 

of " natural selection" based on autocatalytic reactions in sets of prebiotic organic molecules and about the " origin of 
life" beginning with a set of simple organic molecules capable of self-reproduction. The essence of these speculations 
can be summarized as follows: let us consider i > \ different and independent autocatalytic particles, and denote 
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by c\^f < c"^j < ■ ■ ■ < SI] the corresponding bifurcation points. If ca > c^^f , then the stationary density of all 
autocatalytic particles is larger than zero. When the density ca decreases adiabatically to a value lying in the 
interval (c^'t~"^\ c^"'}), j < £, then the autocatalytic particles Xj,Xj+i, . . . ,Xi disappear from the system, and by 



increasing again the density of A particles above cj^^j, there is no possibility to recreate those particles which were 
lost. In this way some form of selection can be realized. When we take into account the stochastic nature of the 
autocatalytic reaction, then we will see that speculations of this kind cannot be accepted. 



B. Forward equation 

Let the random function ^(i) be the number of autocatalytic particles X at the time moment t > 0. The task is to 
determine the probability 



r{m = n\m^Nx}=Pn{t) 



(6) 



of finding exactly n autocatalytic particles X in the system S at time instant t > provided that at t = the number 
of X particles was Nx- Assume that the number of the substrate particles A is kept constant during the whole 
process, we can write for n — 1,2, . . . the equation: 



Pnit + At) = pn{t) 



1 - aNAU At - ^l3'n{n ~ 1) At - At 



i/3'n(7i + 1) +7(n+ 1) 



+aNA{n ~ 1) pn-iit) At + 
and for n — the equation: 

poit + At) ^ Po{t) + jpi{t) At + o{At) 



Pn+i{t) At + o{At), 



where 



kA n/ kx , , 

« = — , " ^ ~^ 7 = fed, 



From these equations it follows immediately that 
dp„{t) 



dt 



{aNA-f3 + j + l3n) n Pn{t) + +aNAin - 1) Pn-i{t) + {(3n + j){n + I) p„+i{t), 



(7) 



(8) 



n= 1,2,..., 



where /3 = (3' /2. If n = 0, then 



Pojt) 
dt 



lPi{i)- 



(9) 



For the sake of completeness we derive the generating function equations 

oo 
n=0 

and 

oo 

5e.p(i,2/) = E{e^«(*)} = ^p„(t) e^y. 

n=0 

By using the equations © and © it is easy to show that g{t, z) satisfies the partial differential equation 

dg{t,z) dg{t,z) d^g(t,z) 
-- -(1 - z) [uNaz - 7) — K y- P(l - z)z 



dt 



dz 



(10) 
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while gexp{t,y) the equation 



= [aNAie^ -!) + (/?- 7)(1 - e — - e . 



(11) 



The initial conditions are g{0,z) — and gexpiO,y) 



— pNxV 



respectively, and in addition it is to note that 



1) = gexp{t,Q) = 1- For many purposes it is convenient to use the logarithm of the exponential generating 
function gexp{t,y)- Therefore, define the function 



K{t,y) = log gexpit,y) 



the derivatives of which at y = 0, i.e. 



dy^ 



Kj{t), J = 1,2, 



are the cumulants of ^ {t) . One can immediately obtain the equation 

- [aNAey - 1) + (/? - 7)(1 - e-)] - /?(1 - e-) 



at 



dy 



dy'^ 



dK{t,y) 



dy 



Now, the initial condition is 

K{0, y) = Nx y, while K{t, 0) = 0. 
In order to simplify the notations define the vector 



(Mt) \ 

Pn{t) 

V : J 



and the matrix 



-Do 


Co 











Bi 


-D, 


Ci 











B2 


-D2 


C2 











B3 


-D3 


C3 



where 



Bn = a {n- 1), 



Dn = n(rt — l + a + 6), 
C„ = (n+ 1) + 

and write the Eqs. and jnj in the following concise form: 

dp(w) 



6 = 



du 



The formal solution of this equation is 



A p(u), where u = (3 t. 



p(u) exp{Au} p(0). 



(12) 



(13) 



(14) 



(15) 



(16) 

(17) 
(18) 

(19) 
(20) 



where the components of p(0) are p„(0) = 5n 



n = 0, 1, ... , and the matrix A is a normal Jacobi matrix. As 



known, the eigenvalues of a normal Jacobi matrix are different and real. If 

> 1/1 > ■ ■ ■ > i/k > ■ ■ ■ 
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are the eigenvalues of A, then the n'th component of the vector p(u) can be written in the form 



p„(w) = ^u;„fce'"=", (21) 

fc=0 



for any n > 0. Since to find the eigenvalues f^, A: = 0, 1, . . . and the coefficients n, fc = 0, 1, . . . is not a simple 

task, we concentrate our efforts only on the determination of stationary solutions of Eqs. © and 



1. Stationary probabilities 

In order to show that the limit relations 

lim Pn{u) = Wn, Vn > 0, n G Zj^ (22) 

u — *oo 

exist, we need a theorem for eigenvalues of the matrix A stating that vq — Q and < for every fc > 1. The proof 
of the theorem can be found in Appendix A. By using this theorem it follows from Eq. H21fl that 

lim p„(u) = w„o, i-e. w„ = Wno, Vn > 0. (23) 

u — *oo 

Since in this case 

dpn{u) w ^ n 

lim ; = 0, V n > 0, 

we can write from immediately the stationary equations 

a [nwr, - {n- l)io„_i] = n[{n + l)w„+i - (n - l)w„] + b [{n + l)w„+i - nw„], (24) 
for all n > 1. Summing up both sides of 124|) from n = 1 to n = fc, we obtain 

a kwk — k{k + l)wk^i + 6 (fc + l)wk+i — b Wi. (25) 
Taking into account the Eq. © we should write that 

1 _ g-kfcl« 



Po{u) = b Pi{u') du' 



oo 



i«o + 5I^ofc e-l-^'^-l", (26) 
fe=i 



and we see that the condition of the stationarity is either wi = 0, or 6 = 0, i.e. 7 = 0. 
If wi = 0, then it follows from 125|) that w„ = for all n = 2, 3, . . ., and consequently, 

00 

EWlk 

Since 

00 00 

Wn =wo = 1, obviously 7 ^ -j — - = 1. 



n=0 fc=l I 

If 6 = 7//3 = 0, i.e. if the autocatalytic particles X do not decay, then wq — and from H25|l one obtains that 

Wk+l = T—r "^k, 1-e. Wn = — |— Wi. 

k + 1 n\ 



Taking into account that X^^^i ~ ^ '^^^ write 



a" 



Wn--0, if 77^0 and wn = -r- if 7 = 0, (27) 

n\ 1 — e"" 



for n = 1, 2, . . ., while for n = we have 

r 1, if 7 7^0, 

wo = <^ (28) 
[0, if 7 = 0. 

For the calculation of factorial moments let us introduce the generating function 

gst{z)^l, if 7 7^0 and gst{z) ^ " \ if 7 = 0. (29) 

1 — e " 
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2. Expectation value and variance 



In order to obtain the equation for the mean value of the number of X particles we use the generating function 
equation Hll|) . Introducing the time parameter u = [3t we have 

— mi(u)[a — b + 1 — mi (u)] — (i [to2 (u) — "t-i (u)] , 

du 

which clearly shows that the kinetic law of the mass action is violated, when the variance V{u) — m2(u) — m\{u) is 
not negligible. For the second moment we get the equation 

P-^ = (a + 5 - 1) mi(u) + [2(a - b) + 3] m2{u) ~ 2 7713(1*), 

du 

in which appears the third moment 7773(77). There are several methods to find approximate solution of ?77i(i), some 
of them were mentioned already in 1]. Here, we do not want to discuss the details, instead we are focussing our 
attention on the properties of the expectation value and variance in stationary state. 

If the decay constant 7 of X particles is zero, then it is easy to show from Eq. H29|l that the stationary value of the 
average number of autocatalytic particles in the system is equal to 



777, = , (30) 



while the second factorial moment is 



,(2) _ " _ „ 



(31) 



consequently, the variance can be written in the form 

V,t = m^^ + m[''^ [1 - m[''^] = t^^^*) [l ~ . (32) 

If 7 ^ 0, then 777^4^ = 777^"*^ = 0, Vfc > 1. 

The first important conclusion drawn from the stochastic model is that the average number of the autocatalytic 
particles in stationary state is different from zero, only when the decay rate constant 7 of X particles is zero. It 
means that there is no "bifurcation" point in the dependence of the average stationary number of X particles on the 
number of A particles, therefore speculations mentioned earlier about the " natural selection" are not supported by 
the stochastic model. 

The second conclusion is connected with the law of the mass action referring to the chemical equilibrium. If 7 = 0, 
then the reversible reaction A + X ^ X + X has to lead to an equilibrium state, in which 

The stochastic model results in an entirely different expression, namely 

(34) 



(st) _ a 



1 - e--^ 
The relative dispersion, i.e. 

Vst 



(st) 



= 1 



clearly shows that the fluctuations of the number of autocatalytic particles become Poisson-like when the number of 
substrate particles is increasing. 
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II. MODIFIED STOCHASTIC MODEL 

The modification of the stochastic model is very simple. We assume that the probability of the reverse reaction 
X + X — > A + X is proportional to the average number of X particles. Therefore, the probability that a reverse 
reaction occurs in the time interval {t, t + At) is nothing else than (3 mi (t) n At + o{At) provided that the number of 
X particles was exactly n at time moment t. Accepting this assumption we can rewrite ^ the equations ((SJ and © 
in the form: 

= -[uNa + + I3mi{t)] n p„{t) + aNA{n - 1) pn^i{t)+ 
+ [7 + /3mi(t)](n + l)p„+i(t), n=l,2,..., (35) 

and 

-^ = 7PiW: (36) 
respectively. One can immediately see that the generating function 

oo 

9exp{t,y) = ^p„(t) e"y 



satisfies the equation 



^hE^ = {NAa{ey-l)-[(3mi{t)+j]{l-e-y)} ^^£|M, (37) 



with gexp{0,y) = e^^^ and gexp{t,0) = 1. In the sequel Nx — 1. Introducing the notations a — Naci/P, b — 7//? 
and u = (3t, it is not surprising that the first moment mi(u) is the solution of the equation 

IIJ- ^ mAu) [a - b - mAu)], 38) 

du 

which is exactly the same as If the initial condition is m„(0) — 1, then 
and one can see that 

a — b, if Na > j/a, i.e. if a > 6, 
0, if Na < j/a, i.e. if a < 6. 
If Na — "f/a — n'^'^\ i.e. if a = 6, then one obtains from H39|l 

TOl(u) = 

1 + u 

The value N^^^'^ is the number of A particles that corresponds to the bifurcation concentration Cbf introduced in the 
rate equation model. 

However, the modified stochastic model takes into account the randomness of the reactions, and hence gives a 
possibility for the determination of the variance of the number of X particles versus time. It can be easily shown that 
the variance of ^ (u) is 



lim TOi(w) = 



D^{^(m)} ~ m2{u) — mflu) ~ mi{t) + mj(u) 



I du' 
2a I —— - 1 

"^lK) 



(40) 



^ In this case we use the notation Pn{t) instead of Pn{t)- 
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In order to prove this relation we need the equation for the second moment m2{u), which can be derived from Eq. 
(|37|l . Introducing the time parameter u = pt we obtain 



dm2{u) 



du 

and this can be rewritten in the form 

dm2{u) 
du 

It is an elementary task to show that 



[a + 6 + mi(w)] mi{u) + 2[a — — mi{u)\ m2(u), 



dmi{u) dlogTOi(u) 
2a mi[u) h 2m2[u) 



du 



du 



m2{u) = mi{u) + 2a mj(u) 



du' 



mi{u') 



and from this we obtain immediately the variance H40|l . 

Taking into account the expression H^W|I for mi(u) the variance of the number of X particles can be written in the 
following form 



D'{u) 



miiu) + mliu) [2 ^(t) - 1] , if Na ^ N^a^\ 



where 



From we can conclude that 



[6u2 + (26+l)u]/(l + u)2, ifNA^N^ 



(bf) 
A ' 



i}{u) 



au a{a — 6 — 1) 



a-b (a -by 



1 - e 



- (a—b)u 



' 0, if Na < N^''^'^ 



A ' 



hm D2{^(u)} = <^ 



b, ifNA^Nl^^\ 
I oo, if iVA > N^A^\ 



(41) 



(42) 
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FIG. 2: Time dependence of the variance of the number of X particles in the case of < Nj^ 



(if) 



In order to have some insight into the nature of the time behavior of the variance of the number of autocatalytic 
particles we have calculated the variance D^{^(w)} versus time curves for different values of the number of A particles 
and for several decay rate constants 7. FIG.|21shows the time dependence of D^{^(w)} when Na < n'^'^'^ at a fixed 
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FIG. 3: Time dependence of the variance of the number of X particles in the case of A'^^ > N at different decay rate constants. 



value of rate constant. We see that the variance versus time curves have a maximum and after that they decrease 
slowly to zero. 

In contrary to this, when Na > N^'^\ then the variance of the number of autocatalytic particles increase mono- 
tonously to infinity. It can be easily shown that 

hm - ' « 



u Na a — 

and so, we can use for large u the asymptotic formula 

D'UH}«^r^ \u, u»l. 

Na a — 

The time dependence of the variance D^{^(u)} in the case when Na > n'^a^ is seen in FIG. 13 for three decay rate 
constants. The variance converges to the infinity linearly with increasing time parameter, if Na > N^'^\ and to a 
constant value b, if Na = N^'^\ so we can say that the fluctuation of the number of X particles in the stationary 
state near the "bifurcation" point alters the possible conclusions based only on the average number of autocatalytic 
particles. 



III. LIFETIME OF THE SYSTEM 

We say that a system is in the state Sj{t), when (,{t) = j. Obviously, the system is live at time instant t > 0, when 
^{t) > 0. Let us define the probability of transition 5; (to) ^ <Sj{tQ + t) by 

P{Sjito+t)\S,{to)} ^p,^j{to + t,to), where t>0. (43) 

Since the process is homogenous in time Pi,j{to + t,to) — pij{t), and using Eq. © we can write that 

dt 

where 



= -(Aj + fij) Pi,j{t) + Xj-i p^j_i(i) + Aij+i Pi.j+i{t), i,j = 0, 1, . . ., (44) 



\,=jNxa, -j[/?(j-l)+7]- (45) 

We see that the first equation is 

^lPi,i{t), (46) 
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and the initial condition is given by pi_j{0) = Sij. The random time 0„ due to the transition iS„(0) — > So{6n), n > 
is the lifetime of the system which has been at i = in the state 5„(0). It is evident that 

v{en >t}^ v{^{t) = o|^(o) = 71} = p„,oW = H„it), 

where Hn{t) is the probabiHty that the hfetime 9^ is not larger than t. The moments of the hfetime are given by 

E{6l^}=/ t'' dH„{t), n>0 and A: = 1,2,... . (47) 
Jo 



A. Extinction probability 

The extinction of a system of state iS„(0) occurs when a transition to the state So{t) is reahzed for any t > 0. We 
define the extinction probabihty by 

hm H„{t) = Lm V 71 = 1,2,..., and since H„{oo) = p„_o(oo) = wq 
in accordance with H28f) we find that 

'1, if 7^0, 

Vn = l,2, ... . (48) 

0, if 7 = 0, 

It is a remarkable result stating that a system being in any of states Snito) at a given time instant to after elapsing 
sufficiently long time t » to will be almost surly annihilated, if 7 ^ 0, and never dies, i.e. the extinction probability 
is zero, if 7 = 0. 

It is to mention that the statement H48|l can be obtained from a nice lemma by Karlin ^ which can be formulated 
in the following way: introducing the notation 

p.=n^, (49) 

J=l J 

where /ij and Aj are non-negative real numbers which are not necessarily equal to (|45|l . it can be stated that if 

n 

K ~ lim pk = 00, then L„ = 1, V n = 1, 2, . . . , (50) 

n — >oo ^ — ^ 

k=l 

while if 

= lim y pfc < 00, then L„ ^ ^ < 1, V7i=l,2,... . (51) 

By using the expressions 145|l we see immediately that 



p.=n 

hence we can conclude that 



/3(j-l)+7 ( P \ r(7//3 + j) 



^.^^ NAa \NAaJ T{-//(3) 



K 



00, if 7 7^ 0, 
0, if 7 = 0, 



and so, applying the lemma we prove the statement H48(l . 
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B. Average lifetime of the system 

To determine the transition probability p„^o(i)j i-G- the probability Hn{t) is not an easy problem. Instead, we show 
how to calculate the average lifetime E{f?„} = t„ . Let us define the parameters 

Si = — , and 6k = hmih^^ if > 1, (52) 
Ml Ml ■ ■ ■ Mfe 

and formulate the following statement called Karlin's theorem Q. // X]^i < oo, then the average lifetime t„ of 
a system containing n autocatalytic particles is given by 

oc n — 1 oo 

k=0 j = l k=j + l 

where pj is defined by i49{ ), and in contrary, if X]fc°=i = oo, then t„ = oo V n > 1. The proof of this statement 
can be found in Appendix B. 

Now, by using this statement we would like to calculate the average lifetime of a system which is in the state 5„ (t) 
at the moment t. Introducing the notations 

and & = J, (54) 

and by using the expressions (|45|l we can write 



where Sj^i is the Kronecker-symbol. Define the sum 



n 1 / ^ 



[j-iy. a^-M_r(& + i) 



L)...(& + j-l) j\ j 7 ^jr(6 + j) 

for n > 1. If 7 ^ 0, then one can see immediately that 

oo 

lim dn — 7 Sj < oo, 

i=i 

i.e. the formula H53|) should be used for the calculation of the average lifetime t„. 
First, determine ti. It follows from Eq. H53I) that 

ri-f:^, = ^^^f:-f^ = - j\{l.b+l-u)du, (57) 

where $(1,6+ 1; u) is the confluent hypergeometric function. The next step is the calculation of the expression 

oo 

Sj{a,b) = Pj X 6k, 

k=] + l 

which can be rewritten into the form: 

j 



i=l k=j+l 



After some elementary algebra we obtain 

.k-3-l 



s,ia,b) = inb + j) Y: T^r 
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and finally we have 



^ 3 = 1 k=j+l ^ ' 

Introducing a new index i = k — j we have 

-J oo 



P^'" ' a + ^ + i)r(6 + j + £ + i)' 

which can be transformed into the expression 

'^^"''^-/3 a^+i J, Z.r(6 + j+£+l) (€+1)! 

By using the identity 



dv. 



r{b + j + i + i) [e + iy. 



Sj (a, b) takes a new form, namely 



1 1 



s,(a,6) = -— ^ / v^-' Ml,b + r,v)-l] dv. 



P a^^ - JO 

Taking into account this formula the expression H58|l can be rewritten in the form 

ji-i 



+ pY.^ v'-'m^,b + j;v)-l] dv, 



which is convenient for numerical calculations. From this equation we see that 

Ti < T2 < ■ • • < T„ < ■ • ■ , 

and we can prove the limit relation 

lim /3 r„ < oo. 



(58) 



(59) 



(60) 
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FIG. 4: The mean value of the system Hfetime vs. 7//3 at A'^^ = 5, 6, 7, provided that the initial state of the system was <Si. 
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It is necessary and sufficient to show that 



(3 hm y Sj{a,b) < oo. 



n — >oo 

J = l 



Since $(1, b + j;v) — 1 is a non-negative, monotonously increasing function of u > we can write immediately that 

/3sj(a,6)< — / dv^ 



ai+^ Jo a j 

and if J > 1, then 

oo ^ oo 

consequently, we obtain the inequality 

n — 1 oo 2 

^ J™o E ^^ («' ^) < E -2 = Y < 

j=i i=i 

and this proves the statement H6Q(I . 

We calculated how the mean value ti depends on the ratio 7//9 at three different values of the number of A particles. 
The results are shown in FIG. 0] We can see that /3ri decreases rapidly with increasing 7//?, and if the values 7//? 
are smaller than 0.4, then we can observe that the larger is the number Na in the system the longer is the average 
lifetime /3ti. 



5 




10 15 
number of initial X particies 



FIG. 5: Dependence of the logarithm of the mean value of the system lifetime on the number of X particles to be found in the system 
at that time moment which the lifetime is counted from. 

That FIG. El shows is rather surprising. The mean value of the system lifetime does not depend practically on the 
number of X particles to be found in the system at that time moment which the lifetime is counted from. We can 
state that the average lifetime of systems controlled by reactions is already determined by several X particles, 
and even a large increase of the number of X particles does not effect significantly on the system lifetime. On this 
basis imagine an " organism" which consists of X particles capable of self- reproduction and self-annihilation. Assume 
that the organism becomes dead if it loses the last X particle. One might think that the greater is the number of X 
particles in the organism the larger is irs average lifetime. Contrary to this conviction, an organism containing, let us 
say 500 particles hardly lives longer than that which contains only 10 particles. By using the values Na — 7, (3 — 0.02 
and 7 = 0.014 we obtain that (3tio ~ 67.21 and Pt^qq w 67.351. It is rather surprising that the increase is only 0.21%. 

In FIG. El we see the dependence of the logarithm of the mean value (3 T20 on the number of A particles at three 
different j/(3. One can observe that the increase of Na results in a rather large lengthening of the average lifetime, 
i.e. the effect of the substrate particles on the process is much stronger than that of the autocatalytic particles. 
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FIG. 6: Dependence of the logarithm of the mean value /3 T20 on the number of A particles at three different values of 7//3. 



IV. CONCLUSIONS 



We assumed the distribution of reacting particles in the system volume to be uniform and introduced the notion of 
the point model of reaction kinetics. In this model the probability of a reaction between two particles per unit time is 
evidently proportional to the product of their actual numbers. By using this assumption we constructed a stochastic 
model for systems controlled by the reactions A + X ^ X + X, X ^ B provided that the number of A particles is 
kept constant by a suitable reservoir, and the end-product particles B do not take part in the reaction. 

We have shown that the stochastic model results in an equation for the expectation value mi{t) of autocatalytic 
particles X which differs strongly from the kinetic rate equation. Further, we found that if the decay constant 7 of 
the particles X is not zero, then the stochastic description, in contrary to the rate equation description, brings about 
only one stationary state with probability 1, and it is the zero-state Sq . It has been also proven that the probability 
of a nonzero stationary state is larger than zero, if and only if the decay rate constant is equal to zero. Consequently, 
the average number of X particles in the stationary state is larger than zero, if only 7 = 0. However, one has to 
underline that this average number is completely different from that which corresponds to the law of the mass action 
of reversible chemical reactions. 

We paid a special attention on the random behavior of the system lifetime, and derived an exact formula for the 
average lifetime. It has been shown that the mean value of the system lifetime does not depend practically on the 
number of X particles to be found in the system at the time instant which the lifetime is counted from. For example, 
the lifetime of a system having 500 X particles at the beginning is larger only by 0.2% than the lifetime of a system 
containing 10 X particles at the time moment t = 0. 



As mentioned already, the eigenvalues of a normal Jacobi matrix are real and different. We would like to prove the 
following theorem: 

Theorem. The v = Q is an eigenvalue of the normal Jacobi matrix A defined by 1^15]) and the other nonzero 
eigenvalues are all negative. 
Proof. Denote by a.y, i, j G Z+ ^ the elements of matrix A. We see immediately that 



i.e. the sum of elements of any column of the matrix A is equal to zero. Strictly speaking, the theorem itself is a 



APPENDIX A: EIGENVALUES OF THE MATRIX A 




(Al) 



^ 2+ is the set of nonnegative integers. 
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straightforward consequence of this property. If v an eigenvalue, then 



U 3 '^'^2 ; 



and xo,xi, . . . ^Xi^ . . . define a nonzero eigenvector 



/ ^0 \ 

Xi 



V • / 



Let y be an arbitrary nonzero vector and form the following expression: 



which can be rewritten in the form 



This equality can be valid only if 



E 



(A2) 



and because y is an arbitrary nonzero vector one can chose its components to be equal to unity. Then taking into 
account the property IjAip one has 



i.e. = is indeed an eigenvalue. 

Now, we show if 7^ 0, then v Since ajj < 0, and > 0, if i 7^ j, it follows from HA2|I that 

As y is an arbitrary nonzero vector, let us chose it so that 



\yj \ — max \yi\ — q > 0, and since y aij = \ajj\, 

i — ^ 



we obtain the inequality 



Taking into account the relation HA3() we have 



which can be valid only ii v < 0. Q.E.D. 



q\a.jj \ > {\a,A + iy)q 



(A3) 
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APPENDIX B: KARLIN'S THEOREM 
Theorem. // J2k=i < oo, then the average lifetime Tn of a system containing n autocatalytic particles is given 

by 

OO 71—1 OO 

fc=0 j = l k=j + l 

where pj is defined by and in contrary, if X]fc°=i ~ then t„ = oo V n > 1. 

Proof. Assume the system to be in the state <S„ at a given time instant t and suppose that the first reaction after 
t occurs at a random time moment t + This reaction can result in a transition to either the state iSn-i or iS„+i 
with probabihties 

and 



An + Pn A,i + p,r, 

respectively. Since t is arbitrary, the equation 

a _ .a , ^" n , -^'i 



is valid with probability 1. Taking into account that 

r{'d„ >t}^l- exp{-(A™ + fin)t} 

one obtains fi'om (IB 21) the recursion 



n+l 



By using the notation 



1 



ifi=l, 



if i > 1. 



(B2) 



Tn-1 + T ^ T„+l, (B3) 



An + /in An + /^n A„ + /i- 

where tq = 0. Introducing the diS'erence w„ = r„ — Tn+i after simple rearrangements we obtain 

C^n - + ^ (B4) 

'^n '^n 

the solution of which can be written in the form: 

_j n— 1 n n 

"..=r+i:f n ^+-n^ 

Z— 1 j—'l-^l J — 1 



and taking into account the identity 
and the relation uj^^ = 'T'i, we have 

(Tn-Tn+l) n-=E'^^-^l- (^5) 

If limn-^oo X]r=i "^i — ^ — OO, then Tn = OO, V n > 1 i.e. the average lifetime of the system is infinite. The proof is 
simple: it is obvious that < Tn+i, therefore, it follows from (|B5|I that '''i ^^r any n, and if n — > oo, than 

Ti must be infinite. Since Tn < t„+i V n > 1, it is evident that t„ = oo, V n > 1, hence the statement is proven. 



If lim„^oo Sr=i = ^ < then one can find a finite real number K such, that nr=i('^*/A'i) < conseque 

lim (t„ - Tn+l) IT ~ = 0' 

and it follows from this that 

oo 
1=1 

Taking into account this relation we can rewrite Eq. ljB5|l into the form: 



Pn 

z=n+l 



^ oo 
{Tu - Tn+l) — ^ 

where pn is defined by H49|) . Introducing the notation 



oo 



2=n 



we have 

Tn+l = Tn + Pn Xn+1, 



the solution of which is nothing else than 



n-l 

Tn = Ti + ^Pj Xj+1, 



3 = i 

and by substituting ri and Xj+i obtain immediately the equation IjBlj). 
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